
% Created 2011-10-27 do 16:51
\documentclass[a4paper, 12pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{fixltx2e}
\usepackage{graphicx}
\usepackage{longtable}
\usepackage{float}
\usepackage{wrapfig}
\usepackage{soul}
\usepackage{textcomp}
%\usepackage{marvosym}
%\usepackage{wasysym}
%\usepackage{latexsym}
%\usepackage{amssymb}
\usepackage{mathtools}
\usepackage{hyperref}
\usepackage{placeins}
\usepackage{pict2e}
\usepackage{subfig}
%\usepackage{cmbright}
\usepackage[dutch]{babel}
\usepackage[a4paper,margin=2.5cm]{geometry}
 \hyphenpenalty=5000
\tolerance=1000
\providecommand{\alert}[1]{\textbf{#1}}
\newfloat{MATLAB code}{h}{}
\usepackage{tikz,pgfplots}
\newcommand{\code}[1]{\textbf{\texttt{#1}}}
\pgfplotsset{compat=newest}
\pgfplotsset{plot coordinates/math parser=false}

\title{Technisch Wetenschappelijke Software: Project Deel 1}
\author{Roel Matthysen \\ s0202264 \\ 1e Master Wiskundige Ingenieurstechnieken}
\date{\today}

\newlength\figureheight
\newlength\figurewidth
\setlength\figureheight{2cm}
\setlength\figurewidth{12cm}


\begin{document}

\maketitle

\section*{Inleiding}
De module geschreven voor deze opdracht rekent de exponenti\"ele functie uit voor een gegeven matrix A. De Taylorreeksontwikkeling van $e^x$ is $\sum_{j=0}^{\infty}\frac{x^j}{j!}$, en voor een matrix A geeft dit
\begin{equation}
e^A=\sum_{j=0}^{\infty}\frac{A^j}{j!}=I+A+\frac{AA}{2}+\frac{AAA}{6}+\cdots
\label{eqn:taylorexp}
\end{equation}
Uitgaande van deze reeksontwikkeling wordt $e^A$ berekend met behulp van de eigenwaardenontbinding van A
\begin{equation}
A=V \Lambda V^{-1}\mbox{\hspace{12pt}met\hspace{12pt}}\Lambda=diag(\lambda_1,\lambda_2,...,\lambda_n).
\label{eqn:eigdecomposition}
\end{equation}
We weten dat $V^{-1}V=I$, en invullen van \ref{eqn:eigdecomposition} in \ref{eqn:taylorexp} geeft dan
\begin{eqnarray*}
e^A&=&I+A+\frac{AA}{2}+\frac{AAA}{6}+\cdots \\
       &=&VV^{-1}+V \Lambda V^{-1}+\frac{V \Lambda V^{-1}V \Lambda V^{-1}}{2} + \frac{V \Lambda V^{-1}V \Lambda V^{-1}V \Lambda V^{-1}}{6}\cdots \\
       &=&V\left(I+\Lambda + \frac{\Lambda \Lambda}{2}+\frac{\Lambda \Lambda \Lambda}{6}\cdots\right)V^{-1} \\
       &=&Ve^{\Lambda}V^{-1}
\end{eqnarray*}
Voorts weten we ook dat $e^{\Lambda}$ gemakkelijk kan worden uitgerekend omdat
\begin{eqnarray*}
e^{\Lambda} &=& I+\Lambda + \frac{\Lambda \Lambda}{2}+\frac{\Lambda \Lambda \Lambda}{6}\cdots \\
		     &=& \begin{bmatrix} 1+\lambda_1+\frac{\lambda_1\lambda_1}{2}+... & & \cdots & 0 \\
		     					 & 1+\lambda_2+\frac{\lambda_2\lambda_2}{2}+... & & \vdots \\
							\vdots & & \ddots & \\
							0 & \cdots & & 1+\lambda_n+\frac{\lambda_n\lambda_n}{2}+... \end{bmatrix} \\
		     &=& diag(\sum_{j=0}^{\infty}\frac{\lambda_1^j}{j!},\sum_{j=0}^{\infty}\frac{\lambda_2^j}{j!},...,\sum_{j=0}^{\infty}\frac{\lambda_n^j}{j!}) \\
		     &=& diag(e^{\lambda_1},e^{\lambda_2},..,{\lambda_n}) 
\end{eqnarray*}

\section{Gevolgde stappen}
In deze sectie wordt een overzicht gegeven van de gevolgde stappen, met bij elke stap een korte uitleg. 

\subsection{Referentieimplementatie van de drie algoritmes}
Vooraleer te beginnen werken in FORTRAN heb ik eerst in MATLAB een referentieimplementatie van de algoritmes \code{taylorexp}, \code{eig} en \code{eigqr} gemaakt, en een testprogramma. Dit was vooral om te kijken of mijn interpretatie van het algoritme klopte, en om te kijken welke LAPACK routines ik juist nodig zou hebben. Wat ik nodig had was iets om eigenwaarden en eigenvectoren te berekenen, een QR-ontbinding uit te rekenen, en een stelsel op te lossen. Deze referentieimplementaties zijn later ook gebruikt om tussenresultaten van mijn algoritmes te vergelijken met de uitkomsten van MATLAB, om fouten te kunnen lokaliseren.
\subsection{Opzoekwerk LAPACK routines, extra argumenten, ...}
De LAPACK routines vermeld in mijn mail waren xGEQRF en xHSEQR. De methode xHSEQR is uiteindelijk niet gebruikt omdat naast de eigenwaarden ook de eigenvectoren nodig waren. Als precisie is telkens gekozen voor re\"ele of complexe dubbele precisie (D en Z).
De gebruikte routines zijn:
\begin{description}
\item[DGEEV] Gebruikt om eigenwaarden en eigenvectoren te berekenen. Deze methode werkte in mijn geval stukken beter dan ZGEEV.
\item[ZGESV] Gebruikt om het stelsel $e^AV=Ve^{\Lambda}$uit het algoritme \code{eig} op te lossen. Vooraleer deze routine toe te passen werden wel transposes genomen om overeen te komen met de conventie van de routine.
\item[ZGEQRF] Gebruikt om de QR-ontbinding van $V$ te berekenen in het algoritme \code{eigqr}.
\item[ZUNGQR] Gebruikt om de matrix $Q$ expliciet op te stellen.
\item[ZTRTRS] Gebruikt om het driehoeksstelsel $RY=Q^{*}$ op te lossen
\end{description}
Om te werken met command-line argumenten heb ik eerst gekeken naar de methodes \code{iargc()} en \code{getarg()}, maar deze werkten enkel met \code{gfortran}. Uiteindelijk heb ik gekozen voor \code{command\_argument\_count()} en \code{get\_command\_argument()}.
\subsection{Structuur van de code schrijven en compileren}
De code is als volgt gestructureerd:
\begin{itemize}
\item Een module \code{methods.f95} die alle algoritmes bevat, plus enkele private subroutines.
\item Een module \code{IO.f95} die het printen en inlezen van matrices regelt.
\item Een wrapper \code{matrixexp.f95} die afhankelijk van de argumenten de verschillende algoritmes aan kan roepen.
\item Een programma \code{test.f95} dat de resultaten voor verschillende testmatrices vergelijkt met de oplossing van MATLAB.
\end{itemize}
Ik ben begonnen met het skelet van de code te schrijven. Dan heb ik de methodes in \code{IO.f95} toegevoegd, en de verwerking van de argumenten in \code{matrixexp.f95}. Ook \code{test.f95} had ik al geschreven voor het implementeren van de algoritmes, zodat elke verandering aan de implementatie onmiddellijk getest kon worden.
\subsection{Schrijven van een testprogramma}
De gebruikte testmatrices zijn voor een deel deze aangereikt in de opgave, en voor een deel zelf aangemaakt met MATLAB. Ik heb eerst het gebruik van system calls overwogen om in het testprogramma de wrapper te kunnen gebruiken. Dit was echter platform-afhankelijk, en werkte ook niet even goed met alle compilers. 

De module \code{test} gebruikt dus de module \code{IO} om de testmatrices in te lezen, gebruikt de module \code{methods} om de algoritmes aan te roepen en drukt de resultaten af op de standaarduitvoer.

Op de resultaten voor de verschillende testmatrices wordt verder ingegaan in sectie \ref{sec:resultaten}.
\subsection{Implementatie van de Taylor benadering}
In het algoritme \code{taylorexp} wordt formule \ref{eqn:taylorexp} expliciet toegepast. De implementatie van deze methode gaf geen noemenswaardige problemen, en de testen gaven goede resultaten.
\subsection{Implementatie van het eenvoudige algoritme}
Bij de implementatie van het algoritme \code{eig} doken al wat meer problemen op, vooral omdat ik ervoor gekozen had om alles meteen in complexe dubbele precisie te implementeren. Een van de grootste was het declareren van dubbele precisie complexe getallen. De functie \code{dcmplx} wordt niet door elke compiler ondersteund. Uiteindelijk was de oplossing om de functie \code{cmplx} te gebruiken met een extra kind-argument.

In de testen gaf het algoritme goede resultaten voor symmetrische matrices, maar hele slechte voor willekeurige matrices. Ik dacht toen dat het aan het verschil in nauwkeurigheid tussen de twee algoritmes lag vermeld in de opgave, dus ik ging verder met de implementatie van het uitgebreide algoritme.
\subsection{Implementatie van het uitgebreide algoritme}
\label{sec:eigqrimplem}
De implementatie van het algoritme \code{eigqr} heeft zeker het langst geduurd, en hier kwam de vergelijking met tussenresultaten van de MATLAB referentieimplementaties zeker van pas. Een moeilijkheid hierbij was dat de eigenvectoren voor complexe eigenwaarden genormaliseerd nog steeds kunnen verschillen in een factor $\pm i$, waardoor de tussenresultaten ook konden verschillen. Een andere moeilijkheid was dat ik eerst de routine ZGEPQF3 gebruikte voor de QR-ontbinding, maar deze gebruikte een kolompivotering, waardoor de matrix $Q$ gepermuteerde kolommen kon hebben.
\newline
\newline Nadat de fouten uit het algoritme waren weggewerkt gaven de testen juiste resultaten voor de voorbeelden uit de opgave, maar enorm grote resultaten voor sommige willekeurige matrices. Door vergelijking met de MATLAB resultaten bleek uiteindelijk dat de methode ZGEEV wel correcte eigenwaarden gaf, maar dat de eigenvectoren soms slechts drie juiste cijfers hadden. Hierna werd ZGEEV vervangen door DGEEV, die wel de eigenvectoren tot op dubbele precisie correct teruggaf. Er moest dan wel een extra private routine toegevoegd worden om de complexe eigenwaarden en eigenvectoren samen te stellen. 
\newline
\newline
Door de berekening van de eigenvectoren te verbeteren was er echter nog geen  verbetering merkbaar in de resultaten voor niet-symmetrische willekeurige matrices. Dit werd na vergelijken met MATLAB tussenresultaten teruggebracht tot onnauwkeurigheden in de QR-ontbinding met ZGEQRF, die zich snel voortplanten als ze daarna in een stelsel terechtkomen. 
Ter illustratie van deze fout een vergelijking van de $R$ uit de QR-ontbinding van $V$, zoals berekend in MATLAB, en zoals berekend met ZGEQRF, voor exact dezelfde $V$ (tot op 15 cijfers nauwkeurig). De $A$-matrix was asymmetrisch met dimensie 4x4.
Uit tabel 1 kan dan worden afgeleid dat de fout voor de drie elementen linksboven quasi nul is, en dat de fout rechtsonder heel snel toeneemt, met uitschieters tot $10^9$. Dit verschijnsel treedt enkel op bij complexe $V$, indien de matrix $A$ symmetrisch is en $V$ dus re\"eel, is de QR-ontbinding wel exact.
\begin{table}[htdp]
\begin{center}
\begin{tabular}{cccc}

      exact  & ${-15.3525} $ &${ -0.8648}$ &   ${-0.8648}$ \\
      0    &  exact &  ${ 4.9486}$ &  ${ 4.9486}$ \\
      0     & 0   & ${4.1903} $ &  ${9.1184} $\\
      0     & 0     & 0 &  ${8.2966} $
\end{tabular}
\caption{$Log_{10}$ van het verschil tussen $R$ berekend met MATLAB en $R$ berekend met ZGEQRF.}
\end{center}
\label{tbl:rfout}
\end{table}%
\newline
Ik zou ook willen vermelden dat zoals ook op het discussieforum gezegd, de uitvoering van ZGEQRF voor complexe $V$ verschilt tussen uitvoeringen, en over het algemeen na de eerste keer uitgevoerd te zijn een kleinere fout geeft.

\subsection{Eenvoudige algoritme revisited}
Door bij de implementatie van \code{eigqr} de eigenwaardenontbinding nader te bekijken, moesten de verbeteringen ook aangebracht worden in \code{eig}. Bij deze methode had de verbeterde nauwkeurigheid van de eigenvectoren wel degelijk resultaat, omdat hierbij geen onnauwkeurige QR-ontbinding moest worden uitgevoerd. De fouten lagen nu in de buurt van de machineprecisie.
\section{Resultaten}
\label{sec:resultaten}
\subsection{Tijdscomplexiteit}
Om de tijdscomplexiteit te testen werd telkens het programma \code{matrixexp} gebruikt, gecompileerd met \code{gfortran}. Aangezien het vooral de bedoeling was om de orde te bepalen, zijn grote willekeurige matrices gebruikt en is er niet zoveel aandacht naar optimalisatie gegaan. De resultaten worden weergegeven in figuur \ref{fig:times}. Uit de figuur blijkt dat de orde van \code{eig} en \code{eigqr} $O(N^2)$ is, en dat de uitvoeringstijden redelijk dicht bij elkaar liggen, met een licht voordeel voor \code{eig} omdat deze de QR-ontbinding niet moet uitrekenen. De curve van \code{taylorexp} heeft duidelijk een andere helling, en de convergentie is eerder $O(N^3)$, waardoor de uitvoeringstijd voor dit algoritme snel toeneemt.
\begin{figure}[htp]
\centering
\includegraphics{img/times.eps}
\caption{De uitvoeringstijd voor de verschillende algoritmes voor oplopende groottes.}
\label{fig:times}
\end{figure}

\subsection{Geheugengebruik}
In Tabel 2 wordt samengevat hoeveel van elke soort variabelen nodig zijn per algoritme, plus een totaal in aantal bytes voor een matrix van grootte NxN. Dit is na een optimalisatiestap, waarbij eerst met behulp van de waarschuwingen van \code{ifort} werd gekeken of er ongebruikte variabelen tussen stonden, en daarna werd gekeken of voor sommige tussenresultaten geen reeds gealloceerd geheugen gebruikt kon worden.

In het eindresultaat is vooral de dominante term van belang, die is voor \code{taylorexp} $24N^2$, voor \code{eig} $136N^2$ en voor \code{eigqr} $152N^2$. De methode met taylorreeksontwikkeling wint dus duidelijk op het vlak van geheugengebruik, en het algoritme \code{eig} doet het iets beter dan \code{eigqr} omdat de extra variabelen voor de QR ontbinding niet moeten worden bijgehouden.
\begin{table}[htdp]
\begin{center}
\begin{tabular}{|rl|c|c|c|}
\hline
 & & Taylorexp & Eig & Eigqr \\
 \hline
 NxN & double complex 	& & 2+1* & 3+1*\\
 NxN & double real		& 3 & 2+3* & 2+3* \\
 Nx1 & double complex	& & 1+1* & 2+1* \\
 10Nx1& double complex & & & 1\\
 Nx1 & double real		& & 2* & 2* \\
 10Nx1 & double real	& & 1* & 1* \\
 Nx1 & integer			& & 1 & \\
 1x1 & integer			& 1 & 3+3* & 3+3* \\
 1x1 & logical			& 1 & 1+1* & 1+1* \\
 \hline
 totaal & bytes			& $24N^2+2$ & $136N^2+132N+8$ & $152N^2+288N+8$ \\
 \hline
\end{tabular}
\caption{Lijst van het geheugengebruik van de verschillende algoritmes. De 10Nx1 vectoren dienen als werkruimte voor de LAPACK routines. De optimale grootte hiervan kan opgevraagd worden door een argument op -1 te zetten, en 10N was voldoende. \\ *=Dit geheugen wordt gebruikt in een subroutine om de eigenwaarden uit te rekenen.}
\end{center}
\label{tbl:rfout}
\end{table}%
\subsection{Nauwkeurigheid}
Om de nauwkeurigheid van de methodes na te gaan werd het testprogramma gebruikt om de verschillende methodes op symmetrische en asymmetrische willekeurige testmatrices van verschillende grootte toe te passen. De resultaten worden weergegeven in figuur \ref{fig:foutensymm}. 
Zoals vermeld in sectie \ref{sec:eigqrimplem} heeft ZGEQRF problemen met de QR-ontbinding van complexe eigenvectoren, dus in het asymmetrische geval stijgt de fout op \code{eigqr} heel snel. De andere twee algoritmes hebben geen problemen met asymmetrische matrices, en het aantal juiste cijfers daalt lineair met het toenemen van de matrixgrootte. Dit lineaire verloop is in het symmetrische geval bij de drie algoritmes zichtbaar.
\begin{figure}[htp]
\centering
\subfloat[Symmetrisch]{
\includegraphics{img/symmetric.eps}\label{fig:symmetric}}
\subfloat[Asymmetrisch]{
\includegraphics{img/asymmetric.eps}\label{fig:assymetric}}
\caption{De maximale fout voor symmetrische en asymmetrische testmatrices van stijgende grootte.}
\label{fig:foutensymm}
\end{figure}
\newline
Het was niet mogelijk om een matrix te construeren waarop \code{eigqr} het merkbaar beter deed dan \code{eig}. De methode die ik gevolgd heb was om eerst een $V$ te maken met conditiegetal $\kappa$ met behulp van een diagonaalmatrix met eigenwaarden en een orthogonale matrix. Het conditiegetal $\kappa$ was dan instelbaar via de diagonaalmatrix. Met deze $V$ werd dan een $A$ geconstrueerd, die dan als input voor het programma werd gebruikt. Door het conditiegetal te verhogen werden de resultaten voor de twee methodes echter allebei snel veel slechter, waardoor het moeilijk was om hieruit conclusies te trekken.
\section{Besluit}
De moeilijkheid van deze opdracht lag niet in het implementeren van de algoritmes, maar in het goed nadenken welke types en welke methodes nodig waren. Het is jammer dat de onnauwkeurigheid van de methodes ZGEEV en ZGEQRF zoveel tijd gekost heeft, maar het zoeken naar de fout was wel zeer leerzaam en heeft zeker meer inzicht verschaft dan een onmiddellijk werkend programma. Het vele debuggen had ook als voordeel dat de geautomatiseerde testen hun doel konden bewijzen. Tijdsbesteding: Ongeveer 40 uur 
\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: